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ABSTRACT 

A numerical scheme is proposed for the solution of the three-dimensional radiative 
transfer equation with variable optical depth. We show that time- dependent ray trac- 
ing is an attractive choice for simulations of astrophysical ionization fronts, particularly 
when one is interested in covering a wide range of optical depths within a 3D clumpy 
medium. Our approach combines the explicit advection of radiation variables with the 
implicit solution of local rate equations given the radiation field at each point. Our 
scheme is well suited to the solution of problems for which line transfer is not impor- 
tant, and could, in principle, be extended to those situations also. This scheme allows 
us to calculate the propagation of supersonic ionization fronts into an inhomogeneous 
medium. The approach can be easily implemented on a single workstation and also 
should be fully parallelizable. 
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1 INTRODUCTION 

Understanding the effect of the radiation field on the ther- 
mal state of interstellar and intergalactic gas is important 
for many areas of astrophysics, and in particular for star and 
galaxy formation. Of special interest for cosmological struc- 
ture formation are the epoch of reionization of the Universe 
(Gnedin & Ostriker 1997, Miralda-Escude et al. 1998), the 
effects of self- shielding on the formation of disk and dwarf 
galaxies (Navano fe 3LeinmeLz 1997, Kepnei el al. 1997) and 



cated structure, and are often of mixed advection-diffusion 
type which makes it very difficult to solve them numerically. 
Besides that, the radiation field in optically thin regions usu- 
ally evolves at the speed of light, yielding an enormous gap of 
many orders of magnitude between the characteristic time- 
scales for a system. 

One way to avoid the latter problem is to solve all equa- 
tions on the fluid-flow time-scale. While there are argum ents 



which seem to preserve causality in such an approach (M 



the abs orption properties of Lya clouda (Katz ct al. 1996, 
Meiksin 1994) . While our knowledge of the physics of large- 
scale structure and galaxy formation has benefited signif- 
icantly from numerical N-body and gas-dynamical models 
(see, e.g., Zhang et al. 1998, and references therein), there is 
very little that has been done to include radiative transfer 
(RT) into these simulations. Challenges seem to abound, not 
least of all, the fact that the intensity of radiation in general 
is a function of seven independent variables (three spatial co- 
ordinates, two angles, frequency and time). While for many 
applications it has been possible to reduce the dimension- 
ality (e.g. to build realistic stellar atmosphere models), the 
clumpy state of the interstellar or intergalactic medium does 
not provide any spatial symmetries. Moreover, coupled equa- 
tions of radiation hydrodynamics (RHD) have very compli- 



halas fc Mihalas 1984), 



even then the numerical sol ution 
is an incredibly difficult challenge (Stone et al. 1992). On 
the other hand, several astrophysical problems allow one 
to follow the system of interest on a radiation propagation 
time-scale, without imposing a prohibitively large number of 
time steps. In the context of cosmological RT, we would like 
to resolve the characteristic distance between the sources 
of reionization. Cold dark matter (CDM) cosmologies pre- 
dict the collapse of the first baryonic objects as early as 
z = 30 — 50, with the typical Jeans mass of order 10 5 M 
(see Haiman & Loeb 1997, and references therein). The cor- 
responding comoving scale of fragmenting clouds is ~ 7kpc. 
If stellar sources reside in primordial globular clusters of 
this mass, then for fib = 0.05 the average separation be- 
tween these objects is Ax ~ 20kpc (comoving). For explicit 
schemes the Courant condition imposes a time-step 
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Thus evolution from z = 20 to z = 3 takes ~ 7.5 x 10 8 /i _1 yrs 
(assuming density parameter flo = 1, and denning the Hub- 
ble constant to be Ho = 100 hkm s _1 Mpc~ ). Substituting 
z — 3 and z = 20 into eq. (Q) yields the total number of time- 
steps in the range 20,000 — 150,000 over the entire course of 
evolution. In other words, to resolve reionization by 10 5 Mq 
stellar clusters, we need to compute ~ few tens of thousands 
of time-steps on average, and for 10 8 Mq clusters the number 
of steps required will be ten times smaller. A 100 3 grid with 
the required resolution will result in computational boxes of 
several Mpc on a side. A full cosmological radiative trans- 
fer simulation with boxes at least this big and for all the 
required timesteps has not been feasible in the past. 

This is far from the only challenge we face. Since any 
two points can affect each other via the radiation field, even 
for a monochromatic problem, we must describe the prop- 
agation of the radiation field anisotropies in the full five- 
dimensional space. Standard steady-state RT solvers, which 
have been widely used in stellar atmosphere models, are not 
efficient in this case. Non-local thermodynamic equilibrium 
(NLTE) steady-state radiative transport relies on obtain- 
ing the numerical solution via an iterative process for the 
whole computational region at once, and is usually effective 
only for very simplified geometries. Any refinement of the 
discretization grid and/or increase in the number of atomic 
rate equations to compute NLTE effects will necessarily re- 
sult in an exponential increase in the number of iterations 
required to achieve the same accuracy. On the other hand, 
the 3D solution of the steady-state transfer equation in the 
absence of any spatial sym metries can often b e obtained 
with Monte Carlo methods (Park & Hong 1998). However, 
these methods demonstrate very slow convergence at higher 
resolutions and are hardly applicable if one is interested in 
following a time-dependent system. 

The change in the degree of ionization in a low-density 
environment occurs on a radiation propagation time-scale 
£r. To track ionization fronts (I-fronts) in this regime, it 
is best to apply a high-resolution shock-capturing scheme 
similar to those originally developed in fluid dynamics. 
One possible approach is the direct numerical solution of 
the monochromatic photon Boltzmann equation in the 5- 
dimensional phase space (Razoumov, in preparation). To 
allow for a trade-off between calculational speed (plus mem- 
ory usage) and accuracy, a more conventional approach is 
to truncate the system of angle-averaged radiation moment 
equations at a fixed moment and to use some closure scheme 
to reconstruct the angle dependence of the intensity at each 
point in 3D space. The method of variable Eddington factors 
first introduced by Auer & Mihalas (1970) has been shown 
to produce very accurate closure for time- dependent prob 



lem s in both 2D (Stone et al. 1992) and 3D ( Umemura et al _ 
1998^~}l owever, to the best of our knowledge, all schemes 
employed so far for calculating the time-dependent variable 
Eddington factor were based on a steady-state reconstruc- 
tion of the radiation field through all of the computational 
region at once, given the thermal state of material and level 
populations at each point. Since advection (or spatial trans- 
port) of moments is still followed on the time-scale of typical 
changes in the ionizational balance of the system, this ap- 
proximation certainly provides physically valid results, as- 
suming that the reconstruction is being performed often 
enough. However, the steady-steady closure relies on the it- 



erative solution of a large system of non-linear equations, 
which becomes an exceedingly difficult problem, from the 
computational point of view, as one moves to higher spatial 
and angular resolution and to the inclusion of more compli- 
cated microphysics. 

The goal of the present paper is to demonstrate that in 
the cosmological context it is possible and practical to solve 
the whole RT problem on a radiation propagation time-scale 
tB, - as opposed to the fluid-flow time-scale - and we present 
a simple technique which gives an accurate solution for the 
angle-dependent intensity in three spatial dimensions. The 
scheme can track discontinuities accurately in 3D and is sta- 
ble up to the Courant number of unity. Since all advection 
of radiation variables is being done at tn, the scheme is 
well tailored to the numerical study of the propagation of 
I-fronts into a non-homogeneous medium with any optical 
depth, and gives very accurate results for scattering. Rather 
than solve the proper chemistry equations applicable to cos- 
mological structures, we adopt a simple toy model described 
below. Similarly, in the current paper we do not make any 
attempt to model the thermal properties of the gas, con- 
centrating just on efficient multidimensional advection tech- 
niques. 

This paper is organized as follows. In Section |2| we 
briefly review the state of numerical RT in the study of 
reionization. We then concentrate on methods for 5D nu- 
merical advection. In Section ^ we describe our numerical 
algorithm and we present the results of numerical tests in 
Section [|. Finally, in Section ^ we discuss the next steps 
towards a realistic 3D RT simulation. 



2 FORMULATION OF THE PROBLEM 

It is believed that light from the first baryonic objects 
at z <; 6 led to a phase-like transition in the ionizational 
state of the Universe. This process of reionization signifi- 
cantly affected the subsequent evolution of structure forma- 
tion (Couchman & Rees 1986) . In detail reionization did not 
happen at a single epoch, with details of 'pre-heating', per- 
colation, helium ionization and other physical processes hav- 
ing been studied in great detail over the last decade (some 
recent contributions include Madau et al. 1997, Haiman & 
Loeb 1997, Gnedin & Ostriker 1997, Shapiro et al. 1998, 
Tajiri & Umemura 1998). 

It now seems clear that the full solution of the prob- 
lem requires a detailed treatment of the effects of RT. To 
complicate matters, by the time of the first star formation, 
the small-scale density inhomogeneities h ad entered the non- 
linear regime (Gnedin & Ostriker 1997), and the medium 



was filled with clumpy structures. The success of cosmolog- 
ical N-body and hydrodynamical models in quantifying the 
growth of these objects (e.g., Zhang et al. 1998) suggests 
that the next step will be to include the effects of global 
energy exchange by radiation. Indeed, there is a need for 
time- dependent 3D RT models as numerical tools for un- 
derstanding the effect of inhomogeneities in the dynamical 
evolution of the interstellar/intergalactic medium. For in- 
stance, the ability of gas to cool down and form structures 
depends crucially on the ionizational state of a whole ar- 
ray of different chemical elements, which in turn directly 
depends on the local energy density of the radiation field. 
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The hydrogen component of the Universe is most likely 
ionized by photons just above the Lyman limit, because (1) 
the cross-section of photoionization drops as v~ A at higher 
frequencies, and (2) the medium will be dominated by softer 
photons, even in the case of quasar reionization (when ion- 
izing photons come mostly from diffuse H n regions) . There- 
fore, we argue that either monochromatic or frequency- 
averaged transfer will be a fairly good approximation in our 
models. 

Recently, the problem of simulating 3D inhomogeneous 
reionization with realistic radiative transfer has attracted 
considerable interest in the scientific community. Umemura 
et al. (1998) calculated reionization from z — 9 to z — 4, solv- 
ing the 3D steady-state RT equation along with the time- 
dependent ionization rate equations for hydrogen and he- 
lium. The radiation field was integrated along spatial dimen- 



sion s us ing the method of short characteristics (Stone et al 
1992^he steady-state solution implies the assumption that 
the radiation field adjusts instantaneously to any changes in 
the ionization profile. One draw-back of this approach, how- 
ever, is in low-density voids where there are probably enough 
Lyman photons to ionize every hydrogen atom, so that the 
velocity of I-fronts is simply equal to the speed of light. Then 
the rate equations still have to be solved on the radiation 
propagation timescale. Besides, implicit techniques in the 
presence of inhomogeneities will become exponentially com- 
plicated, if we want to solve time-dependent rate equations 
for multiple chemical species. 

Norman et al. (1998) and Abel et al. (1998) present 
a scheme for solving the cosmological radiative transfer 
problem by decomposing the total radiation field into two 
parts: highly anisotropic direct ionizing radiation from point 
sources such as quasars and stellar clusters, and the diffuse 
component from recombinations in the photoionized gas. In 
their method the direct ionizing radiation is being attenu- 
ated along a small number of rays, each of which is forced to 
pass through one of the few point sources within the simula- 
tion volume. The diffuse part of the radiation field is found 
with a separate technique which can benefit from the nearly 
isotropic form of this component, for instance, through the 
use of the diffusion approximation. Both solutions are ob- 
tained neglecting the time dependent term in the radiative 
transfer equation, with the default time-step dictated by the 
speed of the atomic processes. 

If reionization by quasars alone is ruled out 
(Madau 1998, however see Haiman & Loeb 1998), then I- 
fronts will be caused by Lyman photons from low-luminosity 
stellar sources at high redshifts. In this case the pressure 
gradient across the ionization zone is more likely to become 
important before the front is slowed down by the finite re- 
combination time. In the present paper we ignore hydro- 
dynamical effects, concentrating on an efficient method to 
track supersonic I-fronts. Our approach is to solve the time- 
dependent RT coupled with an implicit local solver for the 
rate equations. This method gives the correct speed of front 
propagation and it also quickly converges to a steady-state 
solution for equilibrium systems. However, we should note 
that until a detailed comparison is made between explicit 
advection (at the speed of light) and the implicit reconstruc- 
tion (through an elliptic solver), it is difficult to judge which 
approach works best in simulating inhomogeneous reioniza- 
tion in detail. 



Although radiation propagates with the speed of light 
and the intensity of radiation depends on five spatial vari- 
ables, plus frequency and time, the RT equation is inherently 
simpler than the equations of compressible hydrodynamics, 
since its advection part is strictly linear. Non-linearities are 
usually introduced when we are trying to reduce the dimen- 
sionality of the problem. Much of the difficulty thus comes 
from inability to get decent numerical resolution in the 5D 
(or 6D with frequency) space with present-day computers. 

In the current work we have attempted to develop an ef- 
ficient method to describe the anisotropics in the monochro- 
matic radiation field propagating through an inhomogeneous 
medium, which we now describe. 



3 THE NUMERICAL TECHNIQUE 

The classical RT equation (without cosmological terms) 
reads 



+ n ■ V/„ = e v — kJi/, 



(2) 



I dly 

c dt 

where I v is the intensity of radiation in direction n and e„ 
and c„ are the local emissivity and opacity. This equation 
is valid for those problems in which the light-crossing time 
across the computational volume L/c is much smaller than 
the Hubble time, and we can neglect the redshift effects 
within the simulation volume. With proper boundary con- 
ditions one can easily account for the Doppler shift in the 
background radiation. 

The basic idea of our technique is to solve eq. (^) di- 
rectly for the angle-dependent intensity 7(r, n, v, t) at each 
point. The total radiation field at each point is divided into 
the direct (from ionizing sources) and diffuse (due to recom- 
binations in the gas) components: 



E 



E 



_l_ F diB 



(3) 



where i,j, k are the three indices in a rectangular grid. The 
energy density E s J^j k due to direct photons coming from 
point sources of ionization can be easily calculated on the 
3D grid via summation over all sources (assuming that there 
are not too many of them within the volume): 



,i,j,k 



where 



EhuN ph 
Aivcr 2 . , 



i if n,j,k,s < cts 

if ri,j,h,e > cts 



is the optical depth, rij t k,s is the physical distance be- 
tween the current point and the source and t s is the age 
of the source. Note that the rate of emission of photons N p h 
can be modified to allow for variability of sources on short 
timescales. For the diffuse component we use an upwind 
monotonic scheme to propagate ID wavefronts J(r, n, v, i) 
along a large number of rays in 3D at the speed of light. 
Following Stone & Mihalas (1992), we apply an operator 
split explicit-implicit scheme, in which advection of radia- 
tion variables is treated explicitly and the atomic and molec- 
ular rate equations are solved implicitly and separately at 
each point. Unlike Abel et al. (1998), we use rays to track the 
diffuse component of the radiation field and the direct ioniz- 
ing background radiation (streaming into the computational 
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volume). Since we need to draw rays essentially through ev- 
ery grid point in the 3D volume, at first glance this approach 
appears to have very large memory requirements. However, 
efficient placing of the rays can significantly reduce the com- 
putational effort. 



3.1 A uniform and isotropic grid of rays 

At each time-step we are interested in getting a solution 
for the mean radiation energy density and material proper- 
ties on a 3D N 3 rectangular grid. Instead of shooting rays 
though each grid node in 3D, we choose to cover the whole 
computational volume with a separate grid of rays which 
is uniform both in space and in angular directions. Assum- 
ing that the computational volume corresponds to the range 
< x,y, z < 1, we first construct a 2D rectangular base grid 
containing N 2 nodes with coordinates with the origin at the 
centre of the cube 



V3 
1 

2' 



1/2 



N 



1/2 



N 



with 



,N, 



(4) 



where the y3 factor ensures that the entire volume is cov- 
ered with rays, and we shoot rays normal to the base grid 
through all of its grid points. Note that the separation be- 
tween 2D nodes is allowed to be larger than 1/N. To cover 
the whole volume with rays, we then rotate the base grid by 
an angle 8i — tt/2 around the y-axis and by 4>im around the 
z-axis, where the rotation angles are discretized to mimic an 
isotropic distribution of rays (with fewer azimuthal angles 
close to the poles) 



2nm 



N, 



, 1 < m < Na , N& = 2N g cos 6i . 



(•») 



Only those rays which pass through the volume (and 
not all of them do!) are stored in memory. Let the total 
number of all possible orientations of the base grid be 

AUgles ~-Ng(N e -l). (6) 

7T 

The resulting number of rays is significantly smaller than 
TV 3 x iVangios- In fact, with N = 64 and 7V ang i cs = 110 (Ng = 
10), we only require 218,242 rays. 

Now, that we have the grid of rays and the 3D rectan- 
gular mesh, we have to specify the rules of interpolation 
between them. Before we start our simulations, for each 
3D grid point (i,j,k) we also store an array of the closest 
four rays going through its neighbourhood in the direction 
(8i, (j>im)- Since the rays do not pass exactly through 3D grid 
points, we use the values of the intensity on the four closest 
rays in that direction to compute the angular-dependent in- 
tensity. Assume that for the point (i,j,k) the distances to 
the four closest rays going in the direction (6i,(j>tm) are di, 
di, dz and di, respectively. We project the point k) onto 
these rays and read the values of the intensities, which we 



write as 7i, I2, I3 and 1 4. We then calculate the intensity at 
(i,j,k) in the direction {8i,<f>im) according to 



9=1 



dididsdi 



d 



Wqlq 



q=l 



d±d2dsd4 



d a 



(7) 



where the weights w q for q > qo are set to zero if go < 4 rays 
were found in the immediate neighbourhood of k). This 
might be the case close to the edges of the computational 
volume; we shall comment more on this while discussing 
the boundary conditions. The form of eq. (0) was chosen 
specifically because: (1) if a ray labeled q happens to pass 
exactly through the point (i,j,k), then Ii,j,k,i, m = d q ; (2) 
if all four rays encompass the point, min(/ 9 ) < Ii t j t h,t,m < 
max(/ 3 ); and (3) if (i,j,k) happens to be far from all four 
rays, then the resulting intensity will just be an average of 
the four J 9 's. 

At each point on our 3D rectangular mesh we assume a 
piece- wise linear dependence of the intensity I on two angles, 
8 and <t>, 



(8) 



within a spherical rectangular element (or a spherical trian- 
gle adjacent to either of the poles) bounded by the angles 61, 
81-1 in 9 and <f> n and </> n -i in 0, with the rectangular grid 
defined as 



ff f4zi._i , ) > i<;< 
\N e -l 2J' ~ ~ 



2nn 
~Ne~' 



and 



1 < n < N e , 



(9) 



We then integrate the intensity over 4n with appropri- 
ate weights to get the scalar radiation energy density at each 
point (l,n) of the rectangular spherical grid 



Efi? >h = - / i iJik (e, $)dn = -YY iVu^ , 

where the quadrature terms for the integration 



(10) 



1 1-1/2 
n-1/2 



l<8<9 
1 <*<4 



)dSl 



are modified to allow for the non-orthogonal angular grid 
(eq.g. 

Since the advection part on the left-hand side of eq. 
is strictly linear, the simplest way to propagate intensities is 
just to shift wavefronts by one grid zone at each time-step, 
accounting for sources and sinks of radiation. Assuming that 
all discretization points along each ray are strictly equidis- 
tant, the intensity at a point j is updated simply as 



711 - 

n_ ie 



(11) 



Alternatively, one could take special care of the length 
of each ray segment contained within a particular 3D grid 
cell, and use a scheme similar to the third-order-accurate 
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piecewise parabolic advection method (PPA) of Stone & Mi- 
halas (1992). In either case we can track sharp discontinu- 
ities in ID with very little numerical diffusion, and, there- 
fore, our approach is well suited to the calculation of I- fronts. 

3.2 Local chemistry equations 

Since in our calculation all advection of radiation variables 
is performed explicitly, we can solve NLTE rate equations 
separately at each point. This makes it relatively easy to 
implement an implicit solver for all atomic and molecular 
processes. 

To demonstrate the capabilities of explicit advection, in- 
stead of solving the proper chemistry equations for multiple 
species with primordial chemical composition, we have here 
adopted a simple toy model with just photoionization and 
radiative recombination in a pure hydrogen medium. The 
implicit solution of possibly stiff rate equations described 
below can be implemented in a similar manner for more re- 
alistic chemistry models. 

The time evolution of the degree of ionization x c is given 
simply by 

-^2 = (1 - x c )g m E - x c 2 nua, (12) 
at 

with opacity 

k = (1 - x c )nnhpvigm/c. (13) 

Here (Jhi is the photoionization coefficient, E the energy den- 
sity of ionizing radiation and a the recombination coefficient. 
The correct expression for the total emissivity e of the gas 
can be obtained by considering conservation of the thermal 
energy density E t h for matter: 



n+l 



= (l-6)f 1 (x c n )+9f 1 (x c n+1 ) 



AE t - 



E 1 



4neAt = hpvinn Ax c 



(14) 



where all A symbols represent the change of variables during 
one time step. The full recombination coefficient 

a = qi + as (15) 

is the sum of recombination coefficients to the ground state 
(cm) and to all levels above the ground state (ob, the 'case B' 
recombination coefficient), V\ is the frequency just above the 
Lyman limit, and we assume that recombinations in Lyman 
lines occur on a short timescale compared to (leiHCiB)" 1 . 
Similarly, the full emissivity (or gas energy loss through re- 
combinations) is 



e = £i + £B, 



(16) 



where ei/e = cti/a. This simple notation ensures radiation 
energy conservation in eq. (0) for pure scattering of Lyman 
continuum photons (i.e., when ap =0). Eq. (113) does not 
account properly for the number of photons entering the 
volume, so that a large photoionization coefficient gm might 
lead to overproduction of ions. To compensate for this, the 
number of photoionizations inside a 3D grid cell per unit 
time is not allowed to be larger than the number of photons 
actually absorbed inside this cell, i.e. 



< 



E 1 



dx, 

dt ~ hp v\ 
Eq. djjjij) 



"II 



At 



(17) 



are solved separately at each point, given 
the local radiation energy density E. Discretization of 
eq. (|l]) in time yields 
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At 



(18) 



where fi(x c ) is just the right-hand side of eq. ([l2]) and 
1/2 < 6 < 1 for stability. This equation can almost always 
be solved via Newton's method for small enough At. Lin- 
earizing eq. (|lil), we get the (i + l)-th approximation to the 
value of Xr, at time t n+1 : 



n + l,z + l n+l,i . 

Xq — Xc 1" 



dfi(x 



g Xe n+i,i 



1 - 1 



p / n-\-l.i\ 



1 

~ 8At 

n + l, i 



OAt 



(19) 



which can then be iterated. 



3.3 The algorithm 

We start calculations by specifying the initial conditions 
(temperature, degree of ionization, and in the simplest cases 
no radiation field inside the volume) and boundary condi- 
tions (the intensity of radiation entering the volume - all 
outward flux at the edges can freely escape the computa- 
tional box). The inward flux at the boundaries is isotropic 
within 27T and is simply 

where E^ is the average background radiation energy den- 
sity. Since each ray within the volume starts and ends at the 
sides of the volume, we automatically have boundary con- 
ditions for each of the ID advection problems. The density 
field is kept static for all tests in this study, but since the ra- 
diation field is being evolved explicitly at the speed of light, 
one could easily evolve the underlying density distribution 
on a much bigger fluid flow timescale if desired. 

The course of the algorithm at each time step can be 
divided into the following steps: 

• project e„, k v from 3D grid to rays, 

• evolve wavefronts along individual rays by At, 

• project intensities from rays to 3D to reconstruct I v as 
a function of five variables I v {r,Q,cj>), 

• calculate the radiation energy density E(r), 

• solve all chemical equations via an implicit iterative 
technique, separately at each point, and 

• compute new e„, k„. 

At the beginning of each time-step we advect ID inten- 
sities according to eq. ( |ll[ ) along each ray. The numerical 
resolution along each ray is simply set to the resolution 1 /N 
of the 3D rectangular mesh, so that along the ray i the point 
j has a coordinate 



A, 



j ~ 1/2 

N ' 



j = l,...,Ni, Ni = kN, 



where U is the length of the ray segment inside the cube. 
Note that one can have much coarser ID grids along rays, 
speeding up the advection but sacrificing both spatial and 
angular resolution. The advected intensities are then 

projected onto a 3D grid to reconstruct the mean energy 
density at each point using eq. ([7] - [lo|). This operation is 
one of the most demanding from the computational point of 
view, since at each of our N 3 points we have to deal with the 
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angular dependence of the radiation field. When this update 
is done, we solve the matter-radiation interaction equations 
implicitly to compute the local level populations (eq. (Jis|) ) . 
This gives us the 3D distributions of emissivity and opacity 
(eq. ([13] - which are then mapped back to the rays and 
used in the advection scheme at the next time-step. 

This simple scheme which we will refer to as time- 
dependent ray tracing can be used as a stand-alone solver, 
or as a closure scheme for the system of moment equations 
through the use of variable Eddington factors (as in Stone 
et al. 1992). In the absence of any sinks and sources of radia- 
tion, the intensity is conserved exactly along each ray. Since 
the number of rays does not vary with time, our method 
guarantees exact conservation of the radiation energy in 3D. 

Note that if I- fronts do not propagate fast, for instance, 
in the cosmological context where the evolution on the light- 
crossing timescale will normally require many thousand time 
steps, it is possible to update the radiation energy density 
E(r) once every few tens or few hundred time steps, while 
still evolving the intensities and solving all rate equations 
properly at the speed of light. We found that in practice 
this shortcut leads to an increase in speed by a factor of 10 
or even higher without any loss of accuracy. However, it is 
important to compute the energy density properly along the 
edges of I-fronts to guarantee the correct rate of growth of 
ionized regions. 

To further accelerate this computation, we use adaptive 
time stepping to put higher time resolution on the 3D cells 
in the vicinity of I-fronts. Since wavefronts cannot propa- 
gate further than one grid zone during one time step, only 
those cells which just experienced large change in E(r) and 
their immediate adjacent neighbours need a proper update 
of E(r). Depending on the width of I-fronts, typically only 
a few percent of all 3D cells require the new, exact value of 
E(r). 

One advantage of the use of angle- averaged moments 
of radiation is that the advection mechanism is essentially 
reduced to 3D, and it is relatively straightforward to imple- 
ment the multi-dimensional conservation scheme for the lin- 
ear advection part of the moment equations. Then one could 
use a much denser spatial grid for the solution of the mo- 
ment equations, and a relatively course grid for the angular 
reconstruction of the intensity of radiation via ray tracing. 
In practice, however, we have found that the mismatch be- 
tween the spatial resolutions of the moment solver and of 
the ray tracing usually leads to numerical instabilities. In 
what follows, we consider ray tracing only as a stand-alone 
solver. 

Another - perhaps, a better - way of coupling angu- 
lar and spatial variations of the intensity may be an exten- 
s ion of the S pherical Harmonics Discrete Ordinate Method 
( Evans 199£ ). For steady-state transfer problems, instead 
of storing the radiation field, this method keeps track of 
the source function as a spherical harmonic series at each 
point. Although the direct implementation of this technique 
for time-dependent problems is probably not realistic, due 
to the lookback time (i.e. the finite speed of light propaga- 
tion) , the spherical harmonic representation of the radiation 
field might require less storage and might result in smoother 
angular dependence as compared with a pure ray tracing ap- 
proach. 



4 TESTS 

For all of our test runs, except the study of the shadow 
behind a neutral clump in Sec. 4.4, we set up a numeri- 
cal grid with dimensions 64 3 x 10 . The angular resolution 
Aggies ~ 10 2 was chosen to match the equivalent resolution 
of 64 s data points for 5D advection. There are N^ ngleB rays 
passing in the immediate neighbourhood (l/64 3 th of the to- 
tal computational volume) of each 3D grid cell, each ray con- 
taining N p ~ 2/3 x 64 grid nodes. Thus, in 5D we obtain the 
equivalent resolution of N p x 64 3 x ATf ngles ~ 1.1 x 10 9 ~ 64 5 
data points. 

Also, for all simulations we take a 2.5 Mpc (comoving) 
cosmological volume, scaled down to z = 10. All densi- 
ties (from low-density optically thin ambient gas to dense 
clumps) fall in the range Sib = 0.001 — 0.05, assuming a Hub- 
ble constant of 50kms _1 Mpc -1 . This low density contrast 
is chosen to cover the typical range encountered in cosmo- 
logical hydrodynamics, and is ideally suited to demonstrate 
transient features during patchy ionization. The absorption 
coefficients are taken to mimic complete self-shielding (as- 
sumed to be rJjlO) at neutral hydrogen column densities 
higher than 10 18 cm -2 , except where we probe different 
regimes, specifically in Sec. 4.3 where t — 10 corresponds 



and Sec. 4.4 where r = 10 corresponds to 



4.1 An isolated spherical expanding I-front 

A simple problem that would test the ability of the method 
to track I-fronts properly is that of a single, isolated 
Stromgren sphere expanding around a source of ionizing ra- 
diation (see, e.g., Abel et al. 1998). One difficulty of the 
current approach is that rays are drawn in a way to cover 
the whole computational volume uniformly and isotropi- 
cally. For a single point source of radiation this would mean 
that only a small number of rays pass through its neighbour- 
hood. Although intensities along individual rays are strictly 
conserved, there is no guarantee that the energy density has 
the right value far from a source of a specified luminosity. 

At time t = we turn on a point source, which starts 
to blow an expanding Hn bubble around it. Since our algo- 
rithm conserves ID intensities exactly, independently of spa- 
tial or angular resolution, the I-front must propagate at the 
right speed, which can be obtained analytically by equating 
the number of direct ionizing photons to the flux of neutral 
atoms crossing the front. In the absence of radiative recom- 
binations (a = 0), the Hn bubble grows indefinitely with 
the radius 



Mt) = 



et 



Mt<t c 



[3iV ph (t - 2 ^)/47rn H ] 1/3 if t > t c , 



where t c = (A" p h/47rc 3 nH) 1 2 . Parameters used for this cal- 
culation are the diffuse neutral gas density f2b = 0.01 (at 
z = 10) and the central source luminosity N p h = 3 x 10 53 s _1 
(all emitted in photons above the hydrogen Lyman limit). 
The comparison between the numerical speed of the I-front 
and the exact solution from eq. ( |4.l| ) is plotted in Fig. 1. 
The difference between the two always stays within one grid 
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Figure 1. The analytical radius of the spherical Hn bubble 
(solid line) compared to the numerical results (filled data points). 
The error bars indicate the resolution of the grid (zfcAx), while 
the dotted lines give the radius at which the speed of the I-front 
drops below the speed of light. The horizontal dashed line at the 
top shows the radius at which the I-front reaches the boundary 
of the computational volume. 



4.2 An isolated Stromgren sphere in the presence 
of a density gradient 

For this test we put a point source of radiation into a density 
gradient along one of the principal axes of the cube. In the 
absence of diffuse radiation from Hn regions (ai = 0) the 
only ionizing photons come directly from the source in the 
centre, in which case the shape of the ionized bubble would 
be a simple superposition of Stromgren spheres with radii 
Rs (</>) varying wit h the azimut hal angle (f> and given by the 
classical solution (Spitzer 196c) 



So = 47r«B 



n H (r, (j))dr, 



(20) 



where So is the photon production rate of the central source. 
For an exponential density gradient along the y-axis 



logn H = logn H ,i + 



r cos tp + j/o 

Ay 



n-H,2 

«H,1 



(21) 



(riH,i, riH,2 being the hydrogen densities on the opposite 
faces of the cube), the equilibrium Stromgren radius Rs(4>) 
is given by a simple equation 



5*o = 



47TCVB 
~63~ 



[e a+bRsW (b 2 Rl((t>) - 2bR s (<t>) + 2) - 2e a ] , (22) 



where 



lnn H ,i + 



yo ln I rnt ■■ 
na.i 



2 cos < 
Ay 



hi 



Ay 

W H ,1 



and 



We take the physical densities on the opposite faces of 
the volume to be fib.i = 0.001 and Q,b 2 = 0.05, and the 



Figure 3. In order to cover the whole computational volume with 
rays uniformly and isotropically, we set the number of rays inside 
a latitudinal angle interval [9, 9+d9] proportional to cos 6. An odd 
number of rays will necessarily introduce slight left-right asym- 
metry for interpolation between the 3D rectangular grid and the 
radial grid for any non-central cross-section through the volume. 
In this figure we show a 2D schematic representation of the two 
grids employed in our method. To simplify the plot, the source of 
radiation is assumed to be sitting in the centre, and only radial 
rays are drawn here. Due to the finite angular resolution, interpo- 
lation between the 3D mesh and the rays is not symmetric with 
respect to the centre. This shows up in the asymmetry between 
the upper and the lower sides of the ionization contours in Fig. 8. 



luminosity is N p h = 10 5 s . Since we do not solve any real- 
istic atomic rate equations in this paper, we take the value of 
the total hydrogen recombination rate from Hummer (1994) 
for some fiducial temperature (T = 10 4 K). In Fig. 2 we plot 
a time sequence of models, for ionization by a central source, 
with no scattering of Lyman photons (i.e. ai =0). The nu- 
merical solution at t = 4.1 Gyrs appears to be very close 
to the exact one for an equilibrium Stromgren sphere. The 
sharp transition layer between the ionized and the neutral 
regions in the high optical depth regime indicates that, in- 
deed, the scheme introduces very little numerical diffusion 
even when extended to 3D. 



4.3 Ionization in the presence of a UV 
background 

The uniform coverage of the whole volume with rays implies 
that extended sources of radiation will be represented statis- 
tically much better than point sources. A simple test mim- 
icking the evolution of dense clouds in the presence of ion- 
izing radiation is to enclose the computational region in an 
isotropic bath of photons. The simplest way to accomplish 
this is just to set up a uniform, isotropically glowing bound- 
ary at the edges of the cube at t = 0. An effective demon- 
stration of time-dependent ray tracing would be its ability 
to deal with any distribution of state variables within the 
simulation volume. For this test, we set up a density conden- 
sation shaped as the acronym for 'radiation hydrodynamics' 
(RHD), with a density 50 times that of the ambient homoge- 
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Figure 2. The cross-section of the numerical 3D ionization front going through the source of ionization for a model with no scattering, 
shown for six output times. The density of the gas follows an exponential gradient (eq. increasing by a factor of 50 from the upper to 
the lower side of the cube. The dashed line gives the analytic solution (eq. pi) ) for an equilibrium Stromgren sphere without scattering. 
The width of the transition layer between the neutral and fully ionized regions is consistent everywhere with the mean free path of 
ionizing photons. The ionizing source of luminosity N ph = 10 51 s _1 (all emitted in photons just above the hydrogen Lyman limit) is 
marked by the asterisk. The nine contour levels correspond to ionization fractions of x c = 0.1, 0.2, 0.9. 



neous medium. The ambient medium has a constant density 



of fl b 



10" 



and the energy density of the background ra- 



diation is E = 5 x 10~ v\ erg cm Hz sr . Fig. 4 shows 
the result of this run. Most of the low-density environment 
is ionized on the radiation propagation timescale. It takes 
somewhat longer for ionizing photons to penetrate into the 
dense regions. Whether these regions can be ionized on a 
timescale of interest, depends on the ratio of the recombina- 
tion timescale to the flux of background radiation. One can 
easily see ionization 'eating in' to the neutral zone, e.g. in 
the disappearance of the serifs on the letters at late times. 
Note that the width of the ionization fronts does not usually 
exceed one grid zone (Fig. 5). 



4.4 Diffuse radiation from H II regions: shadows 
behind neutral clouds 

Part of the ionizing radiation at high redshifts comes in the 
form of hydrogen Lyman continuum photons from recombi- 
nations in diffuse ionized regions. The following test, simu- 
lating the formation of shadow regions behind dense clouds 
at the resolution 32 3 x 10 2 , was adapted from Canto, Steffen 
& Shapiro (1998). A neutral clump of radius R c = 1.25 Mpc 
(comoving) is being illuminated by a parallel flux F, = 
6.34 x lO'^h' 1 erg cm s Hz 1 of stellar ionizing pho- 
tons (just above 13.6 eV) from one side; here hp is Planck's 
constant. A shadow behind the clump is being photoionized 
by secondary recombination photons from the surrounding 



Hn region (Fig. 6) of physical density tih = 3.7x 10 -2 cm -3 . 
Neglecting hydrodynamical effects, the width Ri of the 
shadow region c an be estimated u sing a simple two-stream 
approximation (Canto et al. 1995): 



l + rf(21n? 



where 



Rc 



and the dimensionless parameter £ is defined as 



4i? c n 2 I a| 



(23) 



(24) 



For £ 2 < 2, recombination Lyman continuum photons 
from the illuminated region will eventually photoionize the 
shadow completely. For £ 2 > 2, radiative losses through low- 
energy cascade recombination photons will stop the I- front, 
forming a neutral cylinder behind the dense clump. Strictly 
speaking, equations (|23[)-(]24|) are valid only for a shadow 
completely photoionized by secondary photons, and should 
be viewed as an approximation to I-fronts driven by scatter- 
ing. 

For this run we modify boundary conditions to include 
recombination photons originating outside the box. Each of 
the rays - starting on any face except the upper side of 
the volume (which goes through the neutral shadow) - car- 
ries the additional intensity of (ei) / (k), where the angular 
brackets denote the average throughout the currently pho- 
toionized gas inside the box. 
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Figure 4. The isosurface of ionization level x e = 0.5 is plotted at six time intervals for the model with an ionizing background 
coming from outside of the cube. The density contrast between the ambient medium and the high-density acronym 'RHD' (radiation 
hydrodynamics) is 50. This simulation demonstrates 5D advection on the radiation propagation time-scale. The numerical resolution is 
64 5 . The effects of shielding are clearly visible during partial ionization. 
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Figure 5. Contour plot of ionization in a cross-section through the data presented in Fig. 4. There are 11 contour lines spaced linearly 
from x e = to x e = 1. Note the sharpness of the ionization boundaries, the shielding, and the fact that the most deeply embedded 
regions take the longest to ionize. 
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Figure 7. The width i?i of the shadow behind a neutral clump 
of radius R c as a function of the dimensionless parameter £ of 
eq. (53). The solid line represents the two-stream analytic solu- 
tion (Canto et al. 1998), assuming an infinitely sharp boundary. 
The points show the results of our 3D numerical model, plotted 
explicitly at x B = 0.5, with the errorbars giving the width of 
the I-front formed by secondary photons. Since the width of the 
neutral core in our numerical model depends on multiple scatter- 
ings in 3D, we would not expect a detailed match; however, the 
agreement between the approximate analytic and the numerical 
solutions is remarkably good. 



In Fig. 7 we plot the radius Ri of the shadow neutral re- 
gion as a function of £ in our 3D numerical models. To play 
with the size of the neutral core, we fix the total recombi- 
nation coefficient but we vary the portion of recombinations 
into the ground state (which would produce more Lyman 
continuum photons capable of ionizing the medium). The 
width of the I-front driven by secondary photons depends 
on the assumed opacity (the optical depth of r = 10 cor- 
responding to the column density of 10 21 cm -2 ). This low 
opacity was chosen to reach equilibrium quicker - equilib- 
rium itself does not depend on the opacity - but it cannot 
be too low otherwise there will be an unnecessary spread of 
the transition zone in the I-front over many grid cells. 

We find a remarkably good agreement between the re- 
sults of our models and the analytic solution, taking into 
account that ionization of the shadow is due to scattering in 
the medium. 



4.5 Diffuse radiation from H n regions: ionization 
of a central void 

To demonstrate the ability of our scheme to handle scatter- 
ing in more complicated situations, we also set up a model 
with ionization of a central low-density void by secondary, 
recombination photons. The void region is surrounded by 
two nested cubes with opposite faces open. The walls of the 
cubes are set to be much denser than the rest of the medium, 
to screen completely the central void from direct ionizing 
photons, and the ionizing UV flux is introduced at all faces 
of the computational volume. The total hydrogen recombi- 
nation coefficient cii + as is again taken for the temperature 



T = 10 4 K. Similar to Sec. 4.4, we vary the amount of scat- 
tering in the medium by changing the fraction q of atoms 
recombining into th e ground state. In reality at T — 10 4 K 
the value q w 0.38 ( Hummer 1994 ) gives a solution in be- 
tween our extreme values of q. 

Similar to the test problem of Section [E^, if qb = 0, 
then the medium will be ionized completely, since there is 
a constant flux of primordial ionizing photons. The speed 
of ionization depends on the values of ati, Qb, gm and na- 
Note, however, that if qi is too high, the I-front will be very 
slow, since a large portion of the original ionizing photons 
are scattered back. On the other hand, if ot\ is too low, the 
I-front will propagate much faster in those regions where 
ionization is driven by primordial photons, but in shadowed 
regions there will be too few recombination photons. Thus, 
it seems that the speed of ionization of the central void will 
be the highest at some intermediate ai. 

In Fig. 8 we demonstrate ionization of the void region 
for models with complete scattering and with no scatter- 
ing at all. The parameters used for this model are E — 
5 x 10 -2C Vi erg cm -2 Hz -1 sr" 1 and the diffuse neutral gas 
density fib = 0.01. 

As expected, for the no-scattering model (q = 0) the 
central region remains neutral, since there is no direct path 
for the ionizing photons. However, for the model which in- 
cludes scattering (q — 1), at least part of the central region 
becomes ionized. This demonstrates that our scheme is per- 
fectly capable of dealing with re-scattering of the ionizing 
photons. Beyond this simple test case, there are many as- 
trophysical situations where progress can be made via nu- 
merical radiative transfer. For example, analytic solutions 
are often used, which are steady-state, and which assume 
a sharp boundary between the neutral and ionized zones. 
Using our numerical techniques it should be possible to fol- 
low general systems with complex density inhomogeneities 
as well as regions of partial ionization. 



5 CONCLUSIONS 

In this paper we have demonstrated that, with existing desk- 
top hardware, it is possible to model cosmological inhomo- 
geneous reionization on a light-crossing time tn, in three spa- 
tial dimensions. Since the photoionization time-scale in the 
low optical depth regime (t < 1) is of order of the light- 
crossing time ta, explicit advection might be a faster method 
in covering at least these regions. Compared to the elliptic- 
type solvers on the fluid-flow time-scale or the time-scale 
of atomic processes, explicit radiative advection produces 
very accurate results without the need to solve a large sys- 
tem of coupled non-linear elliptic equations. The computing 
requirements with explicit advection grow linearly with the 
inclusion of new atomic and molecular rate equations, which 
is certainly not the case for quasi-static solvers (although it 
is feasible that the development of multigrid techniques for 
elliptic equations might actually approach similar scaling). 

Using eq. (|l]) we can see that the entire history of reion- 
ization can be modeled with ~ 10 4 -10 5 time-steps (depend- 
ing on the required resolution), which makes explicit ad- 
vection an attractive choice for these calculations. However, 
the efficiency of the explicit radiative solver has still to be 
explored. Future work should include a detailed compari- 
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Figure 6. Cross-sections of the dense clump (shaded) and the neutral shadow region behind it, for different fractions q = cei/a of 
recombination into ground states. The clump is being illuminated with a parallel flux F* of direct stellar ionizing photons entering 
through one side of the box (marked with arrows), and the shadow region behind it is being photoionized by recombination Lyman 
continuum photons from the surrounding Hn gas. Each model has been evolved to an equilibrium state, correspoding to many passages 
of the wavefront across the computational region. The nine contour levels correspond to ionization fractions of x c = 0.1, 0.2, ...,0.9. This 
test is adapted from Canto et al. (1998). 



son between explicit advection and implicit reconstruction 
(through an elliptic solver), to demonstrate which method 
works best for calculating inhomogeneous reionization. 

As we have demonstrated here, for certain problems, in- 
cluding the propagation of supersonic I-fronts, the Courant 
condition does not seem to impose prohibitively small time- 
steps. In this case the biggest challenge is to accurately de- 
scribe anisotropics in the radiation field, i.e. to solve for 
inhomogeneous advection in the 5D phase space, in the 
presence of non-uniform sources and sinks of radiation. 
Strictly speaking, the storage of one variable at, say, 64 
data points requires about 9 GB of memory, which stretches 
the capabilities of top-end desktop workstations. One at- 
tractive possibility for future exploration is to directly solve 
the monochromatic photon Boltzmann equation in 5D. To 
demonstrate the feasibility of the numerical solution, how- 
ever, among different methods, we here chose to concentrate 
on simple ray tracing at the speed of light. The numerical 
approach we have used is completely conservative and pro- 
duces very little numerical dissipation. 

We want to conclude that despite the high dimension- 
ality of the problem, with a reasonable expenditure of com- 
putational resources, of the type available today, it is possi- 
ble to numerically model many different aspects of the full 
3D radiative transfer problem. Furthermore we feel that the 
methods described here represent a significant and realizable 
step towards the goal of full cosmological RHD. 
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Figure 8. Diagram demonstrating ionization of the void region, which consists of two nested cubes with opposite sides open. The central 
void and the side 'tunnels' have the same density as the surrounding medium, and the density of the walls is set to be much higher to 
completely block all external ionizing photons. In this figure we plot one of the central cross-sections showing contours of ionization at 
five different output times for the model with complete scattering (q = 1), and a quasi-equilibrium configuration for the model with no 
scattering (q = 0) after n = 500 timesteps. The shaded arears represent neutral regions. 
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